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ABSTRACT 

Orthogonal  collocation  using  piecewise  cubic  Hermite  functions  is  used  to 
solve  the  elliptic  partial  differential  equations  arising  from  pseudo¬ 
continuum  models  of  heat  transfer  in  a  packed  bed.  Problems  arising  from  a 
discontinuity  in  the  wall  boundary  condition  and  from  the  semi- inf inite  domain 
of  the  differential  operator  are  discussed.  Comparison  is  made  between  the 
computed  solution  and  experimental  results. 
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SIGNIFICANCE  AND  EXPLANATION 


Models  used  by  engineers  are  often  simplified  in  the  interests  of 
mathematical  convenience,  rather  than  in  accord  with  physical  reality.  A 
typical  case  is  the  omission  of  axial  conduction  from  models  of  heat  transfer 
in  packed  beds,  thus  allowing  an  initial-value  problem  to  be  solved  in  the 
axial  direction,  instead  of  a  boundary-value  problem. 

This  work  demonstrates  the  application  of  orthogonal  collocation,  a 
numerical  method  for  the  solution  of  differential  equations,  to  packed  bed 
models  which  include  axial  conduction.  The  method  is  accurate  and  relatively 
fast,  due  to  the  local  nature  of  the  approximating  functions  used;  these  make 
it  a  good  candidate  for  the  solution  of  "difficult"  problems,  in  which  steep 
gradients  are  involved. 

A  formal  notation  for  the  method  is  given,  which  casts  the  method  in  a 
framework  similar  to  that  of  the  polynomial  orthogonal  collocation  method 
already  familiar  in  the  chemical  engineering  literature  (12,  21).  Some 
tentative  suggestions  for  the  implementation  of  the  method  are  given,  in  the 
form  of  subroutines  supplied  in  an  appendix  to  the  main  text. 
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SOLUTION  OF  PACKED-BED  HE AT -EXCHANGE P  MODELS 
BY  ORTHOGONAL  COLLOCATION  USING 
PIECEWISE  CUBIC  HEPMITE  FUNCTIONS 


Anthony  G.  Dixon 

Introduction 

The  choice  of  a  model  to  describe  heat  transfer  in  packed  beds  is  one  which  has  often 
been  dictated  by  the  requirement  that  the  resulting  model  equations  should  be  relatively 
easy  to  solve  for  the  bed  temperature  profile.  This  consideration  has  led  to  the 
widespread  use  of  the  pseudo-homogeneous  two-dimensional  model,  in  which  the  tubular  bed  i= 
modelled  as  though  it  consisted  of  one  phase  only.  This  phase  is  assumed  to  move  in  pluo- 
flow,  with  superimposed  axial  and  radial  effective  thermal  conductivities,  which  are 
usually  taken  to  be  independent  of  the  axial  and  radial  spatial  coordinates.  In  non- 
adiabatic  beds,  heat  transfer  from  the  wall  is  governed  by  an  apparent  wall  heat  transfer 
coefficient. 

The  earliest  heat-transfer  studies  neglected  the  effective  axial  conduction  term  as 
this  was  expected  to  be  negligible  by  comparison  with  the  bulk-flow  term  in  the  long  beds 
typically  used  in  industry.  Axial  dispersion  was  also  neglected  in  mixino  studies,  and 
experiments  by  Hiby  (1)  confirmed  the  absence  of  axial  back-mixing.  Some  studies  have 
questioned  this  omission  (2,3).  More  recently  it  has  been  shown  (4,5)  that  measurements 
temperature  profiles  in  non-reacting  systems  in  laboratory  packed  bed  heat  exchangers  can 
yield  statistically  meaningful  heat  transfer  parameter  estimates  only  if  the  measurements 
are  made  at  relatively  short  bed  depths,  where  significant  axial  and  radial  temperature 
gradients  are  present.  The  omission  of  axial  conduction  at  such  bed  depths  leads  to 
systematic  errors  in  the  predicted  temperature  profiles,  which  cause  the  model  to  be 
statistically  rejected  when  it  is  fitted  to  data  taken  at  several  bed  depths.  If  the  r- 
is  fitted  depth-by-depth,  the  parameter  estimates  are  found  to  have  a  depth-dependent-, 
noticed  by  De  Wasch  and  Froment  (6).  In  this  case,  they  must  be  regarded  as  lengt1- 
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averaged  values  rather  than  point  values.  Li  and  Finlayson  (7)  argue  that  constant 
asymptotic  values  should  be  used,  as  obtained  from  data  taken  at  long  bed  depths;  this 
would  give  badly-determined  estimates  as  noted  above. 

When  a  chemical  reaction  is  present,  implying  laraer  temperature  gradients.  Young  and 
Finlayson  (8)  have  shown  that  an  effective  axial  dispersion  term  should  be  included,  and 
Mears  (9)  has  given  criteria  for  the  neglect  of  axial  dispersion  which  show  that  increasing 
fluid  velocity  reduces  axial  effects.  This  is  to  be  expected,  since  conduction  through  the 
solid,  a  static  effect,  is  believed  to  be  the  major  contributor  to  axial  effects. 

The  disadvantage  of  including  axial  dispersion  is  that  an  exit  boundary  condition  must 
be  specified,  and  in  cases  where  an  analytical  solution  is  not  available,  a  numerical 
boundary-value  problem  must  be  solved  in  the  axial  direction,  rather  than  an  initial-value 
problem.  This  point  has  been  discussed  recently  by  Jackson  and  co-workers  (10),  who 
present  a  cross-flow  model  of  dispersion  which  allows  downstream  propagation  but  not 
upstream.  However  it  would  not  be  possible  to  use  a  one-phase  model  of  this  form  for  heat 
transfer,  as  axial  back-conduction  cannot  be  represented. 

For  steady-state  heat  transfer  an  elliptic  partial  differential  equation  is  the  result 
of  using  the  one-phase  model.  Previous  studies  (8,11)  have  used  the  orthogonal  collocation 
method  due  to  Villadsen  and  Stewart  (12)  to  determine  the  coefficients  of  trial-f unction 
expansions  in  both  spatial  co-ordinates.  This  method  works  well  when  the  temperature 
gradients  are  moderate  and  few  collocation  points  are  required.  For  steep  profiles, 
however,  such  as  may  be  encountered  at  a  "hot-spot"  in  the  reactor,  many  collocation  points 
may  be  required,  especially  as  the  generation  of  these  points  as  roots  of  polynomials  does 
not  allow  them  to  be  placed  in  the  reaion  of  interest.  Such  a  collocation  scheme  is  a 
olobal  one,  resulting  in  a  collocation  matrix  which  is  large  and  not  usefully  sparse,  so 
that  the  solution  of  the  resulting  algebraic  equations  may  become  costly.  furthermore,  the 
polynomial  collocation  method  has  been  known  to  oscillate  about  the  true  solution  as  the 
degree  of  approximation  is  increased  (13). 


The  answer  to  these  difficulties  lies  in  the  use  of  piecewise  approximants ,  such  as 
cubic  splines,  which  are  in  general  use  in  the  mathematics  literature  (14).  Carey  and 
Finlavson  (15)  have  introduced  a  finite-element  collocation  method  along  these  lines,  which 
uses  polynomial  approximants  on  sub-intervals  of  the  domain,  and  apply  continuity 
conditions  at  the  break -points  to  smooth  the  solution.  It  would  seem  more  straight¬ 
forward,  however,  to  use  piecewise  polynomials  which  do  not  require  explicit  continuity 
equations;  in  this  paper  the  use  of  piecewise  cubic  Hermite  functions  is  considered,  as 
described  by  Prenter  and  Russell  (16). 

The  advantages  of  piecewise  polynomials  are  firstly  that  the  subintervals  may  be 
clustered  in  regions  of  interest,  so  that  an  improved  approximation  may  be  obtained  where 
gradients  are  steep,  and  secondly  that  the  collocation  matrix  is  banded,  allowing  advantage 
to  be  taken  of  this  special  structure,  both  in  work  required  for  decomposition  and  in 
computer  storage  used.  This  second  advantage  was  clearly  demonstrated  in  a  preliminary 
study  on  the  one-phase  model  (17),  where  the  cubic  Hermite  function  method  was  shown  to 
give  an  order-of-magnitude  improvement  in  exit  temperature  profile  over  the  polynomial 
collocation  method  even  when  the  subintervals  were  chosen  to  be  of  equal  length.  A  more 
extensive  invest! oat ion  using  the  one-phase  model  as  a  test  case  is  described  in  the 
present  work. 

The  use  of  two-phase  homogeneous  continuum  models  in  packed  bed  modelling  has  often 
been  avoided  due  to  the  computational  difficulties.  Recently,  Paspek  and  Varma  (18)  have 
found  a  two-n^ase  model  to  be  necessary  to  describe  an  adiabatic  fix* d-bed  reactor,  while 
Dixon  and  Cresswell  (17)  have  shown  that  the  effective  parameters  of  the  one-phase  model 
may  be  interpreted  in  terms  of  the  more  fundamental  parameters  of  a  two-phase  model,  thus 
demonstrating  more  clearly  their  qualitative  dependencies  on  the  operating  and  design 
characteristics  of  the  bed.  '*.Then  two  phases  and  several  species  are  involved,  the 
computational  adi'an^ages  of  the  cubic  Hermite  method  may  be  anticipated  to  be  high. 
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In  this  paper  the  coupled  elliptic  partial  differential  erniations  arising  from  a 


phase  homogeneous  continuum  model  of  heat  transfer  in  a  packed  bed  are  solved,  ar.d 
attempt  is  made  to  discriminate  between  rival  correlations  for  those  parameters  not 
well-established,  by  means  of  a  comparison  with  experimental  results  from  a  previo* 
(4,5). 
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Collocation  using  piecewise  cubic  Hermite  functions 
a)  Boundary-value  problems 

The  interval  [a,b]  is  partitioned  into  n  subintervals  by 

a  =  5.  <  5-  <••■<  5  <5  ,=b.  Then  the  piecewise  cubic  Hermite  functions  are  defined 

12  n  n+1 

for  1  <  i  <  n+ 1  by 


CC)  = 


-2(^i-1)3  3(^i-1>2 

h3  h2 

1-1  1-1 

(  215-5,)3  3(5-5. )2 

'  1  +  r - - - 


(5i-i 4  5  « v 


(5i  4  5  4  w 


otherwise 


(5-5i_1)  (5-5i ) 


(5i_i  <  5  <  V 


i-1 


Vi(5) 


(5-5i)(5i+1-5) 

h2 

1 


(?i  ‘  ?  ‘  W 


otherwise 


The  lenctth  of  subinterval  [5^»5i+1l  is  denoted  by  h^.  The  functions  ^ ,  v 1  are 

restricted  to  the  interval  [5. <5,1  and  $  V  .  are  restricted  to  the  interval 

12  n+1  n+i 

l5n'5n+1!* 

The  piecewise  cubic  Hermite  interpolation  polynomial  to  a  function  u(?)  is 

n+ 1 

u(£>  -  s  <£>  *  I  u(£  )*.  (?)  *  uMC.  • 

n  k  k  k  k 

Hence  s  (£)  interpolates  u  and  its  first  derivative  u*  at  t**e  points  o*  tN» 
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farther ,  let  a . ( £ .  . )  = 

-i  13 


anl  similarly  for  b.  .  and  c.  . . 

-13  -13 


Then 
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The  2n-*-2  unknown  constants  in  3  (£)  are  determined  bv  collocation  at  the  2n  Gaussian 

n 

points  an;  esc  of  the  two  boundary  conditions.  Computational  details  are  given  in  the 
appendix,  and  an  illustrative  example  is  shown  below. 


Example;  comparison  or  so«e  collocation  methods  for  a  boundary-value  problem. 

The  following  problem  a  pood  test  case  (27),  giving  a  "boundary-layer"  type  of 
p  r  o  h  1 « •  -  ; 


u "  (  X )  *■  2 Yxu  *  f  x )  +  2 yu  (  v  )  *  0 


{  1  ) 


v\  j  '  ■ 


.4  . 


u(0 )  =  1 


u(1)  -  e~y 


(la, 1b) 


The  solution  is  u< x)  =  e  T  ,  which  falls  steeply  from  u  =  1  at  x  =  0  at  a  rate 
determined  by  the  parameter  y • 

The  following  methods  are  considered: 

i)  polynomial  collocation,  implemented  using  the  routines  given  in  (21). 
ii)  finite  element  collocation,  implemented  as  described  in  (15). 
iii>  cubic  Hermite  collocation,  as  described  above. 

In  case  (i),  collocation  was  performed  at  the  zeros  of  PN^'^(x),  giving  a  system  of 
N  algebraic  ecruations.  For  methods  (ii)  and  (iii)  the  interval  [0,1]  is  divided  into 
n  subintervals  by  the  set  of  breakpoints  {x^}.  In  the  case  of  finite-element  collocation 
it  was  found  to  be  more  effective  to  collocate  at  only  two  points  per  subinterval,  and  to 
increase  the  number  of  subintervals,  rather  than  to  increase  the  number  of  interior 
collocation  points  in  each  subinterval.  Continuity  equations  must  also  be  included  for 
case  (ii)  at  the  breakpoints?  this  results  in  a  system  of  si2e  3n  +  1,  with  a  bandwidth 


The  equations  for  method  (iii)  in  the  above  notation  are 


(c.  .  +  2yx  ,b,  .  +  2ya.  .).  u.  =  0 
“ID  13-iD  “3-D  -i 

j  *  1,2  ?  i  =  1,2. ..n 


u ( x  )  «  1  ;  u(x  )  =  e 

1  n+1 


These  give  a  system  of  order  2n  +  2,  with  a  bandwidth  of  5.  The  same  breakpoints 
were  used  for  methods  (ii)  and  (iii). 

The  difference  between  approximate  and  exact  solutions  was  examined  at  the  points 
x  "  0,01(0.01)0.10(0.10)0.90;  the  maximum  absolute  error  is  denoted  by  tlell^. 

gone  results  are  shown  in  tables  1  and  2.  The  times  shown  are  relative  and  have  no 
absolj*-*'  meaning.  Both  methods  (ii)  and  (iii)  are  seen  to  be  faster  for  the  same  accuracy 
•  t n  method  (i),  especially  for  the  case  y  =  1500  which  involves  an  extremely  steep 


gradient.  Methods  (ii)  and  (iii)  aave  ♦‘he  san*-  an  expected,  with  method  '  iii) 

heina  faster  due  to  the  fewer  number  of  equations  ♦'eele  }.  ?!v-  performance  of  methods  Mi 

and  (iii)  depended  on  the  breakpoints  use'4,  a  nrr-ur, j f orr  mesv  beinn  *-he  best  ov’oi'~e. 
These  were  chosen  heuristically ,  however  it  should  he  noted  that  s:rply  agtfregat i ng  roart 
near  x  »  0  did  not  always  qive  the  best  results. 
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N  lie  II 
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Time 
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10  0.32 
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i  0(0.05)1.0 
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20  O.S210-3 

4.7 

;  o(o.oi)i.o 

30  0.3910-5 

10.3 

40  0 . 2 1 1 Q -8 

19.4  j  ; 
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0(0.005)0.1 

i! 

(0.01 )0. 3(0.1 ) 

Method  (l) 


il  e  U 


Table  1  y  =  150 


Breakpoint 

set 


Method  (ii)  Method  (iiM 


Time 


Method  (ii)  Method  (iii) 


10 

0.22 

1.7 

20 

0.12 

4.7 

30 

0.13 

10.3 

40 

0.221(1-1 

19.5 

50 

0.59in-3 

32.4 

0(0.001  )0.05  0 . 49  7  p-'J 

(0.01)0.1(0.1) 

1.0 


(0.05)1.0 


e  5  \  =  1500 


The  use  o*  piecewise  bicubic  Hermite  functions  in  collocation  schemes  for  the  solution 
of  elliptic  partial  differential  equations  has  been  described  by  Prenter  (16,20);  a  short 
outline  is  presented  here. 

Partition  the  domain  [a,b]  x  [c,d]  by 


a  =  y,  <  v2  <...<  yn  <  yn+1  =  b 


=  X,  <  X,  <  .  .  .<  X„  <  X„4.,  =  d 


then  the  piecewise  bicubic  Hermite  interpolation  polynomial  to  a  function  S(y,x)  is 


6(y,x)  ~  s  (y ,x) 
n  ,m 


i=i  j=i 


lS(y.,x.)i.(y)?.(x)  +  -5—  ty.,x.)'^.(y);.(x) 
1313  3y  131  3 


2 

l v.  ,x.  )f  ( v) .  (x)  +  (v.,x.)^  (y)0.(x)} 

3x  '  i  3  i  '  3  3y5x  131  3 


This  involves  4(nt1)fm+1)  unknown  constants.  The  Gaussian  points  for  the  x  am!  y 
directions  are  defined  as  before;  combininc;  the  two  points  for  each  of  !vi<Vi+i- 
^xj/kj+1]  gives  four  Gaussian  points  for  the  subrectangle  tVi>vi+1-  *  !xj«xi4.i'- 
Collocation  at  these  points  for  each  subrectangle  yields  a  total  of  4nn  eauations.  1 
should  be  noted  that  at  any  collocation  point  (yi1(,x_..)  only  sixteen  bicubic  croduct 
functions  are  non-zero,  hence  each  collocation  equation  involves  only  sixteen  unknowns . 

Now  let  S  .  =  <9(y.,x.),  ~  (v.,x.),  -ip-  (y  ,x  1,  (v  ,x  ))  and 

-ij  i  3  3y  13  3x  1  3  3y?x  1  3 


— i 3  — i 3 '  Li+l,jf  -i , 3  + 1 '  — i+1,3+1 


Let  a  ( v  ,  x )  =  (  i  .  (  v )  $  .  (  x!  ,v  .  ( v  1  f  (  x)  ,  f  .  ( v  )■;  .  (  x>  .  (  y )  .  (x)>  and 
— 1  i  '  i*i  1  i  1  1  1  ? 


f  r.  \  1  1  or.  linof. 


to  obtain 


(iii)  at  each  corner  both  (i)  and  (ii)  above  may  be  applied,  to  give  four 

conditions.  However  only  three  of  these  will  be  independent,  and  one  must  be 
eliminated,  except  when  there  is  a  corner  discontinuity,  when  an  arbitrary 
decision  must  be  made. 

The  4n  +  4m  +  4  conditions  derived  above  may  be  used  either  to  eliminate  unknowns 
and  thus  reduce  the  size  of  the  system  of  equations  to  be  solved,  as  was  done  in  the 
oriqinal  paper  (16),  or  to  generate  extra  equations.  The  latter  procedure  is  easier  to 
apply  when  mixed  boundary  conditions  are  used,  the  resulting  increase  in  computer  time  used 
hema  offset  by  the  savino  in  programminn  effort. 


c)  Parabolic  partial  differential  equations 

The  solution  of  parabolic  partial  differential  equations  by  the  cubic  Hermite  function 
method  and  by  several  other  methods  has  been  reviewed  by  Hopkins  and  Wait  (23).  They  also 
provide  extensive  program  listings  for  the  solution  of  the  set  of  ODEs  which  results. 


The  packed  bed  heat  exchanger  considered  here  is  that  used  in  recent  experimental 
studies  (3/5)  and  shown  schematically  in  figure  1.  The  long  unheated  calming  section,  (a) 
and  the  heated  test  section,  (b) ,  are  each  considered  to  be  semi-infinite  and  packed  with 
similar  solid  particles.  There  is  a  step  change  in  wall  temperature  at  the  plane  z  =  0. 

The  model  equations  are  well-known,  and  are,  in  dimensionless  form. 


39 


1  39, 


i  lie  _i_  si 

3x  "  Pe,  ,  2  Pe  .  2  +  y  3y 
A  3x  P  3y 


(21 


38 

t—  *  0  at 
3y 


I  3  ) 


8  ♦  0  as  x 


(4) 


8  ♦  1  as  x  +  ® 


(51 


j|  +  (Bi)0  -  (Bi)H(x)  at  y  “1 


(8) 


Here  the  Heaviside  step  function  is  used  to  model  the  change  in  wall  temperature.  The 
downstream  axial  boundary  condition,  equation  5,  was  found  to  be  consistent  with 
experimental  data  in  a  previous  study  (5). 

The  equations  (2)  -  (6)  have  an  easily-determined  analytic  series  solution  in  terms  o 
Bessel  functions  (5);  this  is  therefore  a  good  test  case  upon  which  to  try  out  the 


numerical  method. 


T 


There  are  ♦■wo  points  of  interest  associated  with  the  numerical  solution  of  eauations 
(2)  -  (6):  (i)  the  x-domain  is  infinite  and  (ii)  there  is  a  step-function  in  the  wall 

boundary  condition. 

{ i )  Transformation  of  the  infinite  domain 

There  are  several  wavs  of  dealina  with  an  infinite  domain.  Guertin  et  al  (24) 
chose  perturbation  solutions  of  t*e  model  as  basis  functions.  This  approach  may  be 
difficult  to  extend  to  more  complicated  equations  than  the  non-linear  initial-value 
problems  which  they  considered;  even  so,  the  modeller  must  do  a  considerable  amount  of 
analytical  work  wit*  this  method. 

Birnbaum  and  Lapidus  (2^)  sugaest  the  use  of  polynomials  which  are  orthogonal  over  the 
infinite  domain,  obtaining  these  either  by  using  a  weighting  function  such  as  e  ,  which 
gives  Hermits  polynomials  orthooonal  over  (-00,00),  or  by  transforming  the  infinite  domain 
onto  a  finite  domain,  and  usina  conventional  polynomials  such  as  shifted  Legendre 
polynomials.  The  drawback  to  these  methods  is  that  there  is  no  control  of  the  placement  of 
the  collocation  points,  some  of  which  are  always  included  in  regions  where  the  profile  is 
essentially  flat,  and  are  thus  wasted. 

The  use  of  standard  piecewise  polynomials,  together  with  an  appropriate  transformation 
of  the  infinite  domain,  overcomes  these  difficulties,  provided  some  care  is  taken  in  the 
transformation. 

1  -lx 

Verhoff  and  Fisher  (26)  used  the  form  t  3  —  tan  (— )  in  their  solution  of  the  Graetz 

ji  a 

problem  with  axial  conduction.  Somewhat  neater  equations  result  from  taking 
t  -  tanhfx/i)  giving 
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The  collocation  equations  for  tM«  case  are: 
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Here  is  the  vector  of  unknowns  for  subrectangle  ij. 

The  boundary  conditions  are: 
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The  natural  ordering  of  equations  (If)  nives  a  bandwidth  of  4n  +  15. 

"her  considerin':  t^e  Jeo ree  t'  accuracy  to  he  reauired  from  a  numerical  method,  it  is 
necessary  to  take  into  account  the  potential  uses  of  the  model  equations  being  solved.  It 
would  he  inappropriate  to  recuire  high  accuracv  in  the  present  study,  as  bed  temperature 
rrofiles  are  seldom  accurately  measured,  consecuently  errors  of  1  “(1  in  comparino 
numerics1  md  analytical  results  were  considered  reasonable. 
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f +  ,pi''  - ?i  p<t>  (ii> 

Initial  lv  it  nay  appear  simplest  to  let  1=1,  and  use  t  =  tanh  x.  However  this  results 
ir  '■'re  recipe  x  >  3  in  *he  bed  being  rannpti  orto  a  small  interval  [0.905,1].  This  is 
undesirable  for  two  reasons: 

1>  If  subintervals  are  reouired  in  the  bed  downstream  of  x  -  3,  it  becomes 

difficult  to  place  the  breakpoints  in  the  t-domain 

2^  if  axial  oradier.ts  are  present  downstream  of  x  =  3,  they  will  become  very  sharp 

in  the  t-domain,  since  —*  B  - t—  so  4“  *  00  as  t  ♦  1. 

5t  n-t2,  3*  3t 

It  can  be  seen  that  for  a  fixed  choice  of  {t^},  different  choices  of  a  will  locate  the 
induced  partition  ;x.'-  in  different  physical  marts  of  the  bed.  Thus  a  must  be  chosen 
so  that  the  collocation  points  in  the  t-dom  in  are  placed  in  a  way  that  makes  physical 
sense  in  the  x- domain.  The  appropriate  value  is  found  by  trial  and  error;  an  empirical 
rule  suaoested  by  experience  is  to  take  a  *  Pe^. 

The  effect  of  an  inappropriate  value  for  a  is  shown  in  figure  2.  It  should  be  noted 
that  although  Verhoff  and  Fisher  used  a  =  1  throughout,  their  computations  were  made  for 
sufficiently  short  beds  to  ensure  good  results;  laraer  values  of  a  would  be  reauired  for 
temperature  profiles  further  downstream,  which  would  be  of  interest  in  the  presence  of 
reaction . 

( i i )  Step-change  in  wall  boundary  condition 
It  was  anticipated  that  the  discontinuity  in  wall  temperature  at  x  =  0,  and  the 
resulting  steep  local  gradients,  would  lead  to  a  locally  poor  approximation  which  might 
•  ,|>-».  r*  V*err-;e  efforts  further  downstream.  it  was  soon  found  that  mesh  refinement  in  t^e 
-'vi.-i1  cirecMon  1nr"pfi  t*e  results  »-onsi  derably  ov«r  the  vise  of  an  eoually-spaced  mesh, 
v-«rrta-  -.ee*’r  refinement  in  the  radial  direction  had  little  effect,  and  a  fairlv  coarse 
-  »"  -*  imI  7*;r-v.  wa»-  al  wavs  *rund  to  adeauate. 


The  mesh  refinement  was  carried  out  by  heuristically  searching  for  good  breakpoint 
distributions,  the  final  choice  (corresponding  to  m  =  10)  being 

{t. }  =  {-0.1,  -0.05,  -0.025,  0,  0.025,  0.05,  0.1,  0.2,  0.6}  . 

Ary  further  refinement  led  to  improved  values  only  at  extremely  short  bed  depths. 

The  effect  of  the  wall  temperature  discontinuity  may  also  be  mitigated  by  an 
alternative  implementation  of  the  wall  boundary  condition  to  that  described  in  the  previous 
section.  The  two  equations  at  each  interior  boundary  node  (1,tj),  j  =  2,3, ...,m,  are 
dropped,  together  with  one  equation  at  each  of  (1,-1)  and  (1,1).  These  2m  eauations 
are  replaced  by  application  of  the  boundary  conditions  at  two  points  within  each 
subinterval  on  the  line  y  »  1.  The  Gaussian  points  are  not  necessarily  the  optimal 
choice,  but  were  used  in  the  absence  of  any  other  guideline. 

For  the  case  of  a  uniform  mesh,  the  Gaussian  point  implementation  was  an  order-of- 
magnitude  improvement  over  the  breakpoint  implementation.  When  the  refined  mesh  was  used, 
the  two  methods  gave  essentially  the  same  results,  except  at  low  bed  depths  near  the  wall, 
where  the  Gaussian  point  method  was  slightly  better.  Consequently  it  was  the  method  used 
in  the  rest  of  the  work. 

The  reason  for  the  above  differences  is  not  clear  but  appears  to  lie  in  the  wall 
temperature  specification.  No  difference  between  the  methods  was  found  when  applied  to  the 
centre-line  condition,  and  reduction  of  Bi  to  lessen  effects  of  the  discontinuity  also 
greatly  reduced  the  advantage  of  the  Gaussian  point  method.  It  may  be  that  the  boundary- 
condition  is  more  effective  when  "spread  out"  using  two  points  instead  of  being  reinforced 
by  applying  two  conditions  at  one  point. 


The  comparisons  were  made  at  y  -  0.1,0. 2. ..1.0  for  each  of  5  bed  depths:  x  = 
1.33,2.61,4.00,5.33,6.67.  The  parameter  ranges  covered  were 

Pe,  :  0.25  -  20.0 

A 

Pep  :  2.5  -  120.0 

Bi  :  0.5  -  8.0  . 

The  refined  mesh  given  above  was  used  in  the  axial  direction;  the  radial  mesh  used  n 
=  3,  {y. }  =  {0.25,0.5,0.75}.  It  was  found  possible  to  use  the  same  mesh  throughout  in  the 
t-donain,  due  to  the  freedom  to  vary  the  scaling  factor  a.  This  avoided  heuristic 

searching  for  a  new  mesh  for  each  new  set  of  parameters.  Automatic  mesh  generation  was  not 

felt  to  he  worthwhile,  in  view  of  the  extra  costs  involved  and  the  relatively 
underdeveloped  state  of  the  art  (27). 

The  1°C  criterion  was  met  in  all  cases  except  for  Peft  =  20  using  this  method.  For 

that  case  the  centre-line  discrepancy  rose  to  approximately  3“C  at  lower  bed  depths. 

Presumably  this  error  could  be  eliminated  by  taking  higher  order  approximations.  Typical 
computation  times  to  produce  one  set  of  solutions  (i.e.,  5  bed  depths)  were  approximately 


1.5  sec  on  a  UNIVAC  1110  computer,  a  reasonable  compromise  between  cost  and  accuracy. 


Two-phase  continuum  model 

Two-phase  continuum  models,  in  which  the  solid  particles  and  their  associated  stagnant 
fillets  of  fluid  are  regarded  as  a  continuous  pseudo-solid  phase,  are  to  be  preferred  to 
the  more  traditional  cell  model  of  a  particle  bathed  in  fluid,  which  does  not  allow 
conduction  from  particle  to  particle.  In  previous  studies,  such  models  have  been 
simplified  by  considering  a  one-dimensional  model  only  (28).  This  study  considers  the  full 
equations,  which  are,  in  dimensionless  form: 
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A  semi-analytical  solution  to  these  equations  was  derived  bv  Dixon  and  Cresswell  (19),  wv*-> 
then  matched  the  fluid  phase  temperature  profile  to  the  one-phase  model  profile  to  obtain 
explicit  relations  between  the  parameters  of  the  two  models. 

The  numerical  solution  to  the  system  of  equations  (12)  -  (17)  parallels  that  of  the 
one-phase  model  almost  exactly,  with  longer  computation  times  due  to  the  increased  size  of 
the  collocation  matrix  and  its  bandwidth.  Typical  computation  times  to  produce  fluid 


solid  temperature  profiles  at  each  of  five  bed-depths  were  3-4  seconds. 


As  the  fluid  was  air,  Pr  =  0.72;  the  bed  voidage  i  was  taken  as  a. 4.  T-e  ratios 

krs/kg  and  kp/kg  were  related  using  the  formula  of  Zehner  and  Sch  Kinder  (30'. 

Preliminary  sensitivity  tests  showed  that  the  parameter  h’u,,,.  had  very  little  e'fert  on 

the  profiles;  in  the  presence  of  reaction  this  parameter  would  he  more  important. 

It  is  clear  that  by  adjusting  the  parameters  of  this  model  by  nonlinear  regression 

excellent  fit  to  the  experimental  data  could  be  obtained.  The  value  of  such  a  proce  iur 

rather  dubious,  however,  and  it  is  more  useful  to  use  the  model  to  obtain  Qualitative 

information  about  the  Quantities  Pe  P»  (•»>  and  Pi  ,  which  are  poorlv  deter-i-e 

rf  af  s 

in  the  literature. 


--1- 


a )  Pe  ,{*) 

This  quantity  is  almost  certainly  a  weak  function  of  the  tube-to-part icle 
diameter  ratio  {  although  the  exact  dependence  is  not  well-established. 

Examination  of  reported  literature  data,  especially  the  results  of  Kunii  et  a! 

(31)  at  extremely  hicih  Ke,  suggests  taking  Pe  ^(°°)  -  8.0. 

s 

b)  Pe  _(*) 

af 

Recent  work  of  Hsiang  and  Haynes  (32)  shows  the  commonly-taken  value  of 
Pe^ff00)  =  2.0  to  be  questionable  in  the  (dt/dp)  range  considered  here.  It  was 
decided  to  use  t^e  relationship  Peaf  21  Pea'  which  results  from  the  model 
matching  of  (19)  when  Re  >  20-30.  The  values  of  Pea  were  poorly  determined 
from  measurements  made  at  bed  exit  alone  (5),  so  estimates  from  data  which 
included  profiles  at  z  =  0  were  used  (4).  These  values  led  to  great 
improvements  in  the  slopes  of  the  predicted  profiles. 

c)  Bi 
_ s 

Values  of  Bis  may  lie  between  the  theoretical  lower  bound  derived  by  Olbrich 
(29): 

Pi  >  (2.12)  (~) 

S  d 

P 

and  31^  *  correspond! no  to  no  thermal  resistance  between  solid  and  wall. 

If  equation  (20)  is  used  to  predict  Bi„,  f-en  it  is  necessary  to  tahe  Bi,.  -  10(ln 
for  a  aood  fit,  as  shown  in  'iqure  3  for  a  typical  case.  However,  it  should  be  noted  that 
equation  (20)  underestimates  the  values  o'  Bi,(»  Pi  I  found  i -  (5).  This  is  probably  due 
to  the  unreliable  correlation  used  for  Hu A.c,  as  pointed  out  in  (4).  If  the  experimental 
estimates  of  Bi  are  used  instead  o'  equation  (20),  then  values  o*  P.i  in  the  ranoe  in  - 

21  are  needed.  Some  o'7  these  results  are  shown  in  fiuures  -  7. 
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This  result  indicates  that  only  precise  determination  of  Bi^  will  allow  any 
conclusions  to  be  drawn  on  Bis#  since  these  parameters  may  be  mutually  varied  over  fairlv 
large  ranges  and  similar  results  obtained 4 

The  computed  temperature  profiles  in  figures  4-7  show  good  general  agreement  wit1- 
experiment;  some  deviation  is  apparent  in  the  centre  of  the  bed  for  Re  «  73  and  224. 
Agreement  is  improved  if  lower  values  of  Pe^l®)  are  used,  but  there  is  no  justification 
for  this  in  the  literature. 

The  semi-analytical  solution  given  in  (19)  may  also  be  compared  to  experimental 
results  using  the  same  parameter  results  as  for  the  numerical  method.  This  is  shown  in 
figure  8  for  the  case  Re  *  430.  Similar  results  to  the  numerical  method  are  obtained  at 
longer  bed  depths;  at  short  bed  depths  the  approximation  is  poor,  due  mainly  to  the 
inadequacy  of  the  one-point  radial  collocation  method  used. 
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Conclusions 

The  orthogonal  collocation  method  usina  piecewise  cubic  Wermite  polynomials  has  been 
shown  to  give  reasonably  accurate  solutions  at  low  computing  cost  to  the  elliptic  partial 
differential  equations  resulting  from  the  inclusion  of  axial  conduction  in  models  of  heat 
transfer  in  packed  beds.  The  method  promises  to  be  effective  in  solvina  the  nonlinear 
equations  arising  when  chemical  reactions  are  considered,  because  it  allows  collocation 
points  to  be  concentrated  where  they  are  most  effective. 

The  *luid-phase  temperatures  predicted  from  a  two-phase  pseudo-homogeneous  model  were 
shown  to  give  reasonable  aareement  with  experimental  measurements,  without  explicitly 
adjusting  the  model  parameters.  It  was  demonstrated  that  more  refined  experimental 
measurements  will  be  needed  to  determine  the  parameters  of  the  model;  in  particular,  the 
solid  and  fluid  phase  wall  Biot  numbers  were  mutually  adjustable. 
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NOMENCLATURE 


specific  interfacial  surface  area  (m*1) 
fluid  specific  heat  (kJ/kg°C) 
pellet  diameter  (m) 
tube  diameter  (n) 

2 

superficial  mass  flow  rate  (kg/m  s) 

o  O 

apparent  interphase  heat  transfer  coefficient  (w/m  C) 

2° 

apparent  wall  heat  transfer  coefficient  (w/m  C) 

2° 

wall-fluid  heat  transfer  coefficient  (w/m  C) 

2° 

wall-solid  heat  transfer  coefficient  (w/m  C) 

2° 

fluid-solid  heat  transfer  coefficient  (w/nr  C) 
axial  effective  conductivity  (w/m°C) 
radial  effective  conductivity  (w/m°C) 
axial  conductivity  of  fluid  phase  (w/m°C) 

axial  conductivity  of  solid  phase  (w/m°C) 

radial  conductivity  of  fluid  phase  (w/m°C) 

radial  conductivity  of  solid  phase  (w/m°C) 

molecular  conductivity  of  fluid  (w/m°C) 
pellet  conductivity  (w/m°C) 
length  of  packed  test  section  (m) 
tube  radius  (m) 
radial  coordinate  (m) 
bed  temperature#  one-phase  model  ( °C ) 
fluid  phase  temperature  ( °C ) 

solid  phase  temperature  (°C) 

temperature  of  calming  section  wall  (°C) 
temperature  of  test  section  wall  (°C) 
superficial  fluid  velocity  (m/s) 
axial  co-ordinate  (m) 


Dimensionless  parameters 


Bi 

apparent  wall  Riot  number,  hwR/k 

Bif 

fluid-wall 

Piot  number, 

hwfRArf 

solid-wall 

Biot  number, 

hwsR/krs 

h . 

subinterval 

length 

m  number  of  axial  subintervals 


-31- 


number  of  radial  subintervals 

interphase  heat  transfer  group,  aR2h/krf 

interphase  heat  transfer  group,  aR  h/krs 

apparent  wall  Nusselt  number,  hwdf/kg 

fluid-solid  Nusselt  number,  hfsdp/kg 

fluid-wall  Nusselt  number,  hwfdpAg 

effective  axial  Peclet  number,  Gcpdp/ka 

effective  axial  Peclet  number  (based  on  R),  GcpR/ka 

effective  radial  Peclet  number,  GCpdp/kr 

effective  radial  Peclet  number  (based  on  R),  GcpR/kr 

axial  fluid  Peclet  number,  GCp^p/^af 

axial  fluid  Peclet  number  (based  on  R) ,  GcpR/kaj 

radial  fluid  Peclet  number,  Gcp^p/'crf 

radial  fluid  Peclet  number  (based  on  R),  GCpR/krj 

asymptotic  value  of  perf  as  Re  ♦  “ 

asymptotic  value  of  Peaj  as  Re  »  “ 

Prandtl  number,  wc  /k 
P  g 

Reynolds  number,  Gd^/u 
normalised  radial  co-ordinate  (r/R) 
transformed  axial  co-ordinate 

dimensionless  fluid  temperature  (Tbf“T0 )/(Tw“To ' 
dimensionless  solid  temperature  (Tbs-T0 )/(Tw-T0 ) 
normalised  axial  co-ordinate  (z/R) 


a  axial  scaling  factor 

e  bed  voidage 

$  cubic  Hermite  basis  function 

ip  cubic  Hermite  basis  function 

i 

0  dimensionless  bed  temperature,  *Tb_T0  >/(Tw-Tg  * 

y  viscosity  of  fluid  (kg/ms> 

p  density  of  fluid  (kg/m^) 

£  a  general  independent  variable 
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Appendix :  computations 1  detail* 

For  collocation  at  twe  Gaussian  points  wit!' in  a  subinterval  #  the  vectors  a^.-,  hj  ^ , 
c ^  ierer.d  onlv  on  and  are  oiven  by 


.1  2 

, .  i , 

h. 

l 

1  2 

,  i 

in  =  !I  +  IT?' 

(1  +  /?’ 

1  2 ' 

2  '  IT?' 

7?  " 

ii2  "  !I  '  IT?' 

1 1  ”  Tf* 

h 

a 

12  ' 

1  +  2 

I  +  IT?' 

(-1  - 

7? 


h 

II1 


’’  4’ 


b 

-il 


— i2 


1  1 

h.  '  IT?' 

l 


IT?1 


i  i 

“ '  IT?'  “ '  I7?! 
1  1 


£n  =  '■ 


2/T  ~/3  -  1 
.  2  '  h.  ‘ 


2/3 
.  2  ' 


1  -  /I, 


,2/3  /I 


2/3  /I  +  1 


— i2 


h. 

i 


These  vectors  are  supplied  by  a  FORTRAN  subroutine: 

CALL  BASWTS  (IP,  H,  V)  . 

This  subroutine  computes  the  weiohts  for  u*ID*(xj,)  and  u^ID' (x.^  * !  i-e.  ID  is  an 
indicator-variable  set  to  0,1  or  2  for  the  function  or  its’  first  or  second  derivative 
respectively.  The  leorith  of  the  subinterval  containinq  x^ ,  and  x^2  Is  supplied  to  the 
subroutine  in  H,  and  the  weiohts  are  output  in  the  2x4  arrav  V  in  the  form  e.q. 

%1> 

K2> 

for  17  ~  0,  ar*f1  similarly  for  IT  =  1,2. 


Tb*»  u**' 


disrretise  boundary  value  problem*  result*  in  a  system  of 


A£  =  F(_u ) 

where  _u  is  the  vector  of  unknowns.  The  matrix  A  corresponds  to  the  linear  part  of  th 
differential  operator  (and  may  be  SO).  If  the  entire  problem  is  linear  then  F (_u ) 
reduces  to  a  constant,  b,  and  solution  is  straight-forward;  otherwise  iterative  methods 
must  be  used.  Computational  advantage  may  be  taken  of  the  fact  that  A  is  a  band  matrix 
if  the  equations  are  ordered  naturally,  some  thought  may  have  to  be  given  to  the  placemen 
of  the  boundary  equations  to  obtain  minimum  bandwidth. 

The  band  matrix  A,  of  order  N  with  Ml  subdiagonal  elements  and  total  bandwidth 
MT,  is  stored  as  an  N  x  MT  array.  The  matrix  is  factorised  in  the  form  A  =  LU,  the 
lower  triangle  being  stored  in  array  XM  and  the  upper  overwriting  A,  by  the  subroutin 
BANDET;  subroutine  BANSOL  is  then  called  to  solve  the  banded  system  LUu  =  J}.  The  routin 
are  simplications  of  those  given  in  (22),  where  further  details  on  the  storage  method  are 
given. 

Once  the  vector  _u  has  been  obtained,  the  solution  at  arbitrary  X  is  found  using 
the  formula 


u  ( x)  =  (a^x)  •  ui> 


where  x^  <  x  <  x^+^. 

The  vector  a^(x)  is  given  by  the  subroutine 

CALL  INTWTS (ID,X,XI,XJ , VECT )  with  ID  =  0  ; 

x  is  supplied  in  X,  x^  in  XI  and  x^+^  in  XJ.  Then  on  output  VECT  is  the  4- 
vector  a^tx). 

If  is  required,  INTWTS  is  called  with  ID  =  1;  ID  =  2  gives  c^(x). 


-3A- 


Note  B ASWTS  is  a  special  form  of  INTWTS  where  advantage  is  taken  of  the  need  to  obtai 


weights  only  at  t'ne  Gaussian  points  of  each  subinterval. 

It  is  easy  to  see  that  the  elements  in  a^j(y,x),  a^ ( j+1 (y ,x) ,  _a^+1  j(y,x)  and 
£^+1 1  j+i  (y  ,x)  may  be  generated  from  those  in  a_^(y)  and  a^(x).  Hence  _A^  ^  <  y  *  3< )  may 
also  be  generated  from  _a^(y)  and  _£j(x)  and  the  product  o  is  defined  Ly 


^jiy.x)  =  a^(y)  oa^(x) 


Similar  relations  hold  for  B^ly.x)  etc.  More  generally  a  16-element  vector  VC  may  be 
generated  from  any  two  4-element  vectors  VA  and  VB  via  the  o-product,  and  a  subroutine  to 
do  this  is:  CALL  OPROD (VA,VB, VC) . 

Collocation  within  subrectangle  ij  requires  A  ^  for  k,i  =  1,2  given  _a^  and 

a  for  k,L  =  1,2.  Subroutine  PRDWTS  calls  OPROD  four  times  to  achieve  this:  a  typical 
2*  2 

o  9 

calling  sequence  is  e.g.  to  obtain  weights  for  — -  on  subrectangle  ij 

3y 


CALL  B ASWTS  (2, HI, VI) 

CALL  B ASWTS  (0,HJ,VJ) 

CALL  PRDWTS  (VI,VJ,VYY) 


Here  VI  and  VJ  are  2x4  matrices  as  described  above  and  VYY  is  a  4  x  16  matrix 

(cyy  )  ’ 

-ij”  | 

(cyy  ) 

-i  j21  i 

r~yy  J 

-ig12 


<cyy  > 

' — i  j  2  2 
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«■  Y  J 


*  'll.' 


IF 


SUPROl’TINE  PANPFTt  NP, N,MTP,MT,M1p,v  1  ,  A,  XM ,  IP,  IV,  IFAIL ) 
IMPLICIT  rxMJPLP  PRECISION  (A-H,0-71 
DIMENSION  IP(NP)  ,  ACID, NIP)  ,X“(NP,M1P) 

IFAIL  =  0 
L  =  Ml 

DO  30  I  =  1,M1 
IJ  =  Ml  +  2  -  I 
DO  25  J  -  IJ,MT 
25  A(I,J-L)  =  A ( I , J ) 

L  =  I.  -  1 
JL  =  NT  -  L 
DO  30  J  =  JL,MT 
30  A(  I ,  J)  =  0 
L  =  Ml 

DO  100  K  =  1,N 
X  =  At  K ,  1 ) 

I  -  K 

IF  (L.LT.N)  T,  =  L  +  1 
IF  (K+1 .GT.L)  GO  TO  40 
IK  »  K  +  1 
DC  35  J  =  IK,L 

IF  (DABS(  A(  J,  1  )  )  .LE.DAPStX  )  )  GO  TO  35 
X  =  A(  J,  1  ) 

I  =  J 

35  CONTINUE 
40  IP f K )  =  I 

IF  (X.NE.0)  GO  TO  50 
IFAIL  -  1 
IV  =  K 
RETURN 

50  IF  (I .EQ.K)  GO  TO  70 
DO  65  J  «  1 ,MT 
X  -  A(K,J) 

A(K,J)  =  A(I,J) 

65  A(  I ,  J)  =  X 
70  IF  (K+1. GT.L)  GO  TO  100 
IK  =  K+1 
DO  00  I  =  IK,L 
XM(K,  I-K)  =  At  I,  1  )/A(K,  1  ) 

X  =  XM(X,I-K) 

DO  SO  J  «  2 , MT 

RO  A{  I ,  J-1  )  =  At  I ,  J)  -  X*A(F,J) 

DO  A(I,MT)  -  0 
100  CONTINUE 
PETUPN 
r.JI; 


>•  v  «X  " 


pi'skctine  '< anppl ( nt , n , mtp ,  mt  ,  m  i  p  ,  m  i  , r.xm.ip.p) 

IMPLICIT  POC'-LE  PRECISION  (A-H,0-7) 

DIMENSION  A(VP,M-"P),  XMrvP,"1p),IP(«T)/P(Vnl 
pp  40  K  =  1  ,N 
I  =  IP(K) 

IF  (I.EO.K)  oC  TP  IP 
X  «  3(K) 

P(K)  »  Pd1 
fo  *■  x 

10  IF  (L.LT.N)  I  «  L  +  1 
▼  p  (X+1.GT.L'  G^  TO  4P 

IX  =  K+1 

PO  2P  I  =  IK,L 

X  =  XM ( K , I -K ) 

20  B(I)  =  Ed)  -  X*B(K) 

4  0  CONTINUE 


no  ii  =  i,n 
I  =  N+1-H 
X  =  B(I  ) 

IX  =  1-1 

IF  (L.EO. 1 )  GO  TO  65 
DC  60  K  -  2,L 
60  X  =  X-A(I ,K)*P(IW+K) 
65  P(I)  =  X/A(I,  1  ) 

IF  (L.t.T.MT)  L  =  L+1 
on  CONTINUE 
PETVPN 


fukpouti: ne  raswts/io,h,v> 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-7) 
DIMENSION  V(2,4) 

S3  =  rSQRT ( 3  - 0P0 ) 

IF  (IP-1)  1,2.3 

1  V(1,1)=0. 5+2. 0/(3. 0*S3> 

V(1 ,3)=1.0-V(1 , 1 ) 
V(1,2)*(1.0+1.0/B3)/12*H 
V ( 1 ,4 )=<1 .n/S3-1.0)/12*H 
GO  TO  4 

7  v(  i ,  i  )«-i  .n/H 

V(1,3)=-V(1,1) 

V( 1 , 2 )=P . 5/S3 
V(1,4)  =  -Vf  1,2) 

GO  TO  4 

3  V< 1 , 1 )  =  -2.n*S3/(H**2  ! 

V ( 1 , 3 )=-V (1,1) 

V(1 ,2 )«-(1 .0+S31/H 
V(1,4)='1."-E3)/H 

4  .:-(- l)**(IP+D 
V ( 2 ,2 )=J*V (1,4) 

VI 2,4 1,2) 

V<  2,3)  =-J+V (1,1) 

V  f  2 , 1 )  =  -J+V( 1,3) 


‘“'VV'-l-An r.- 


SUBROUTINE  INTWTS ( ID ,X , XI , X J, VECT) 

IMPLICIT  DOUBLE  PRECISION  (A-H.O-Z) 

H-XJ-XI 
Z-(X-XI)/H 
Y=(X-XJ)/H 
IF  (ID-1)  1,2,3 

1  VECT(3)-(3-2*Z)*Z**2 
VECT(1 )»1-VECT(3) 

VECT ( 4 )  -H  *Y*  Z*  *2 
VECT(2)-H*Z*Y**2 

GO  TO  5 

2  VECT( 1 )-6/H*(Z**2-Z) 

VECT{3)— VECTd) 

VECT(4)-Z**2+2*Z*Y 
VECT(2)=2*Z«Y+Y**2 
GO  TO  5 

3  VECT(1 )=6/H**2*(2*Z-1 ) 

VECT ( 3 J--VECT  ( 1  ) 

VECT(2)-(2*Z+4*Y)/H 

VECT(4)-(4*Z+2*Y)/H 

5  RETURN 
END 

SUBROUTINE  PRDWTS(V1 ,V2,V) 

IMPLICIT  DOUBLE  PRECISION  (A-H.O-Z) 

DIMENSION  VI  (2,4)  ,V2(2  ,4)  ,V{4, 16)  ,VA(4  )  ,VB(4)  ,W(  16  ) 
IPT-0 

DO  10  J=1,2 
DO  10  1-1,2 
IPT-IPT+1 
DO  15  K-1 ,4 
VA(K)-V1(I,K) 

15  VB(K)»V2(J,K) 

CALL  OPROD  ( V A ,  VB ,  W  ) 

DO  10  K-1, 16 
10  V(IPT,K)=W(K) 

RETURN 

END 

SUBROUTINE  OPROD( VA,VB , V ) 

IMPLICIT  DOUBLE  PRECISON  (A-H.O-Z) 

DIMENSION  VA( 4 ) ,VB( 4 ) ,V( 16 ) 

JJ-1 

DO  10  J-1,2 
DO  10  1-1,2 
K-2*(  1-1  )  +  1 
L-2* (1-1 )+1 
V( JJ)«VA(K)*VP(L) 

V(JJ+1 )-VA(K+1 )*VB(L) 

V(JJ+2)-VA(K)*VB(L+1  ) 

V(  JJ+3)=VA(K+1  )*VB(L+1  ) 

10  JJ-JJ+4 
RETURN 
END 
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